package scu.maqiang.matlabcases;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarCFunc;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Complex;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.IterSCSolver;
import scu.maqiang.numeric.SCMatrix;
import scu.maqiang.numeric.SRMatrix;

public class Scattering {

	public static void main(String[] args) {
		String meshFile = "matlabMeshes\\Scattering.msh";
		Mesh2T3 mesh2D = new Mesh2T3().readMesh(meshFile, MeshFormat.msh);
		mesh2D.toTecplot("matlabMeshes\\ScaterringMesh.dat");
		
		FES2T31 fs = new FES2T31(mesh2D);
		SRMatrix RA = new SRMatrix(fs.GetNdof());
		// RC = new SCMatrix(fs.GetNdof());
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, RA);
		ScalarFunc coefC = (xy, label, param) -> -3600.0;
		fs.assembleMass(new ScalarFunc[] {coefC}, null, BVPType.CONSIST_MASS, RA);
		SCMatrix LHS = new SCMatrix(RA);
		
		double k = 60.0;
		ScalarCFunc fluxFunc = (xy, label, param) -> new Complex(0.0, -k);
		fs.assembleBoundaryMass(new ScalarCFunc[] {fluxFunc}, null, BVPType.CONSIST_MASS, LHS, 5, 6, 7, 8);
		//LHS.sort();
		//System.out.println(LHS.getElement(77, 76));
		//System.out.println(LHS.getElement(78, 78));
		
		fs.applyBC_MBN(LHS, Direct.All, 1, 2, 3, 4);
		ScalarCFunc fixFunc = (xy, label, param) -> new Complex(-Math.cos(k * xy[0]), Math.sin(k * xy[0]));
		Complex[] RHS = Complex.zeros(fs.GetNdof());
		fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{fixFunc}, null, 1, 2, 3, 4);
		//System.out.println(LHS.getElement(45, 45));
		//System.out.println(Complex.toString(RHS));
		IterSCSolver solver = new IterSCSolver(LHS);
//
		Complex[] x = Complex.zeros(fs.GetNdof());
		solver.PCGSSSOR_2(RHS, x, 1.5, 1);
		//for(int i = 0; i < fs.GetNdof(); i++) {
		//	System.out.println(x[i].image);
		//}
		mesh2D.toTecplot("matlabResults\\Scattering.dat", x);
	}

}
